This is an Rmarkdown notebook. If you’re viewing it in a web browser, you can click on Code->Download Rmd in the top right corner, then open it in RStudio and to run/edit the code. If you don’t use RStudio you can also just copy the code into your editor of choice.
We’ll be using a variety of tidyverse packages as well as sf and leaflet. If you don’t have these packages installed, uncomment the following code and install them now. (It’ll take a little while.)
# install.packages("tidyverse")
# install.packages("sf")
# install.packages("leaflet")
#load just the tidyverse for now
library(tidyverse)
Let’s create some very simple point data.
five_points <- data.frame(var1 = 1:5, var2 = as.integer(c(3,5,4,4,1)))
five_points
Call ggplot() to make a canvas
Define aesthetics with aes() - these map your data to spatial properties like the x and y dimensions, or visual aesthetics like color, fill, alpha etc
Add a geom - here, geom_point()
five_points %>%
ggplot(aes(x=var1, y=var2)) +
geom_point()
%% is the modulo function)five_points <- five_points %>%
mutate(var1_odd_or_even = ifelse(var1%%2==0,"even","odd"),
var2_odd_or_even = ifelse(var2%%2==0,"even","odd"))
five_points
five_points %>%
ggplot(aes(x=var1, y=var2,color = var1_odd_or_even)) +
geom_point(size = 5)
Color is for lines (including outlines of shapes)
Fill is for areas (filling in shapes)
but our points are filled in based on the color aesthetic?!
the default point shape is shape 19 (“circle”) - just a solid circle with the color taken from the color aesthetic; shape 21 (“filled circle”) is a circle with a fill and an outline, which can be specified separately
Generally, you can either map the aesthetic to a data point or specify its value outside of the aes() function. If you do neither, a default will be used.
Aesthetic mapping/specification can be done in the original ggplot() call, or within individual geoms
Good general reference on specifying aesthetics:
https://ggplot2.tidyverse.org/articles/ggplot2-specs.html#point-1
five_points %>%
ggplot(aes(x=var1,
y=var2,
color = var1_odd_or_even,
fill = var2_odd_or_even)) +
geom_point(size = 5,
shape = 21,
stroke = 2)
Longitude and Latitude are (on a 2d map) x and y coordinates
Let’s get some coordinate data from the NYC facilities database
https://data.cityofnewyork.us/City-Government/Facilities-Database/ji82-xba5
public_libraries <- jsonlite::fromJSON("https://data.cityofnewyork.us/resource/ji82-xba5.json?factype=PUBLIC%20LIBRARY")
public_libraries %>% str()
'data.frame': 221 obs. of 34 variables:
$ uid : chr "015fec883329f3481f1711a39386a5ee" "05375983dd54899915412bc855795f75" "0616b9a31fe2db14793561377be01c0e" "08d71649b166891363a946ebb0b3b029" ...
$ facname : chr "MORRIS PARK LIBRARY" "NORTH FOREST PARK" "58TH STREET LIBRARY" "FLUSHING" ...
$ addressnum: chr "985" "98-27" "127" "41-17" ...
$ streetname: chr "MORRIS PARK AVENUE" "METROPOLITAN AVENUE" "EAST 58 STREET" "MAIN STREET" ...
$ address : chr "985 MORRIS PARK AVENUE" "98-27 METROPOLITAN AVENUE" "127 EAST 58 STREET" "41-17 MAIN STREET" ...
$ city : chr "BRONX" "FOREST HILLS" "NEW YORK" "FLUSHING" ...
$ boro : chr "BRONX" "QUEENS" "MANHATTAN" "QUEENS" ...
$ borocode : chr "2" "4" "1" "4" ...
$ zipcode : chr "10462" "11375" "10022" "11355" ...
$ latitude : chr "40.84806105240" "40.71118635740" "40.76223213860" "40.75778232920" ...
$ longitude : chr "-73.85697658820" "-73.85367190670" "-73.96932374730" "-73.82887651450" ...
$ xcoord : chr "1023819.55661000000" "1024817.14797000000" "992747.99653700000" "1031658.10613000000" ...
$ ycoord : chr "248281.25916200000" "198414.74758100000" "216979.94069500000" "215403.54794900000" ...
$ bin : chr "2045398" "4076687" "1037165" "4114282" ...
$ bbl : chr "2041260040.00000000000" "4032070022.00000000000" "1013130005.00000000000" "4050430011.00000000000" ...
$ commboard : chr "211" "406" "105" "407" ...
$ council : chr "13" "29" "4" "20" ...
$ censtract : chr "25200" "72900" "11203" "85300" ...
$ nta : chr "BX37" "QN17" "MN19" "QN22" ...
$ facgroup : chr "LIBRARIES" "LIBRARIES" "LIBRARIES" "LIBRARIES" ...
$ facsubgrp : chr "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" "PUBLIC LIBRARIES" ...
$ factype : chr "PUBLIC LIBRARY" "PUBLIC LIBRARY" "PUBLIC LIBRARY" "PUBLIC LIBRARY" ...
$ capacity : chr "0" "0" "0" "0" ...
$ optype : chr "Public" "Public" "Public" "Public" ...
$ opname : chr "New York Public Library" "Queens Public Library" "New York Public Library" "Queens Public Library" ...
$ opabbrev : chr "NYPL" "Public" "NYPL" "Public" ...
$ overlevel : chr "City" "City" "City" "City" ...
$ overagency: chr "New York Public Library" "Queens Public Library" "New York Public Library" "Queens Public Library" ...
$ overabbrev: chr "NYPL" "QPL" "NYPL" "QPL" ...
$ datasource: chr "nypl_libraries" "qpl_libraries" "nypl_libraries" "qpl_libraries" ...
$ facdomain : chr "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" "LIBRARIES AND CULTURAL PROGRAMS" ...
$ schooldist: chr "11" "28" "2" "25" ...
$ policeprct: chr "49" "112" "18" "109" ...
$ servarea : chr "Local" "Local" "Local" "Local" ...
We need latitude and longitude to be numeric, not strings
public_libraries <- public_libraries %>%
mutate(latitude=as.numeric(latitude),longitude=as.numeric(longitude))
public_libraries %>%
select(facname,latitude,longitude) %>%
str()
'data.frame': 221 obs. of 3 variables:
$ facname : chr "MORRIS PARK LIBRARY" "NORTH FOREST PARK" "58TH STREET LIBRARY" "FLUSHING" ...
$ latitude : num 40.8 40.7 40.8 40.8 40.7 ...
$ longitude: num -73.9 -73.9 -74 -73.8 -73.7 ...
Map the longitude and latitude to the x and y aesthetics
public_libraries %>%
ggplot(aes(x=longitude,
y=latitude)) +
geom_point()
public_libraries %>%
select(latitude,longitude) %>%
summary()
latitude longitude
Min. : 0.00 Min. :-74.24
1st Qu.:40.66 1st Qu.:-73.98
Median :40.72 Median :-73.92
Mean :40.35 Mean :-73.25
3rd Qu.:40.77 3rd Qu.:-73.85
Max. :40.90 Max. : 0.00
public_libraries %>%
filter(latitude<39|longitude>74) %>%
select(latitude,longitude,facname,addressnum,streetname,address,city,boro)
public_libraries <- public_libraries %>%
filter(!latitude==0,
!longitude==0)
public_libraries %>% nrow()
[1] 219
public_libraries %>%
ggplot(aes(x=longitude,
y=latitude)) +
geom_point()
Make it a bit prettier
Color the points and rename the legend with scale_color_manual()
Get rid of axes and background with theme_void()
#define some system colors - a named vector of colors will map the colors to the values of the color aesthetic (backticks let you use non-syntactic R names)
#color values can be names from R's built-in colors or hex codes
library_system_colors <- c(`New York Public Library`="red3",
`Brooklyn Public Library` = "orange2",
`Queens Public Library` = "purple3")
public_libraries %>%
ggplot(aes(x=longitude,
y=latitude,
color = overagency)) +
geom_point() +
scale_color_manual(values = library_system_colors,name = "System") +
theme_void() +
coord_equal()
sf (simple features) packagehttps://r-spatial.github.io/sf/
library(sf)
Simple features is a set of standards for defining two-dimensional geometries, used by various GIS systems, building up from x-y coordinates within a coordinate reference system (crs).
POINT
LINESTRING
POLYGON
MULTIPOINT
MULTIPOLYGON
sf package in R lets you:
store geometry data as a list-column in a data frame or tibble, alongside other data columns
import data into the R sf format from various GIS formats
manipulate sf geometries (concatenate, subtract, find overlaps, etc)
You can make an sf object directly from a data frame:
five_points_sf <- five_points %>% st_as_sf(coords = c("var1","var2"))
five_points_sf %>% print()
Simple feature collection with 5 features and 2 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 1 ymin: 1 xmax: 5 ymax: 5
CRS: NA
var1_odd_or_even var2_odd_or_even geometry
1 odd odd POINT (1 3)
2 even odd POINT (2 5)
3 odd even POINT (3 4)
4 even even POINT (4 4)
5 odd odd POINT (5 1)
There is a special geom for sf objects: geom_sf
When you add geom_sf to a plot (and the sf package is loaded), ggplot will plot the geometry column - plotting method depends on the sf geometry type (point, polygon, etc)
five_points_sf %>%
ggplot() +
geom_sf(
aes(color = var1_odd_or_even,
fill = var2_odd_or_even),
size = 5,
shape = 21,
stroke = 2
)
But for mapping we usually import geometry data from an external dataset.
Formats:
GEOJSON
ESRI shapefile
Many others (st_drivers(what = "vector"))
NYC geography resources:
https://www1.nyc.gov/site/planning/data-maps/open-data/census-download-metadata.page
#read geojson format directly from download link
pumas_2010_geojson <- st_read("https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/NYC_Public_Use_Microdata_Areas_PUMAs_2010/FeatureServer/0/query?where=1=1&outFields=*&outSR=4326&f=pgeojson") %>%
st_as_sf()
Reading layer `OGRGeoJSON' from data source
`https://services5.arcgis.com/GfwWNkhOj9bNBqoJ/arcgis/rest/services/NYC_Public_Use_Microdata_Areas_PUMAs_2010/FeatureServer/0/query?where=1=1&outFields=*&outSR=4326&f=pgeojson'
using driver `GeoJSON'
Simple feature collection with 55 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS: WGS 84
pumas_2010_geojson %>% print()
Simple feature collection with 55 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID PUMA Shape__Area Shape__Length geometry
1 1 3701 97928415 53226.65 MULTIPOLYGON (((-73.89641 4...
2 2 3702 188993596 106167.48 MULTIPOLYGON (((-73.86477 4...
3 3 3703 267645154 305260.19 MULTIPOLYGON (((-73.78833 4...
4 4 3704 106217077 47970.15 MULTIPOLYGON (((-73.84793 4...
5 5 3705 122483734 68697.43 MULTIPOLYGON (((-73.87046 4...
6 6 3706 43887042 51825.93 MULTIPOLYGON (((-73.87773 4...
7 7 3707 42281072 37374.60 MULTIPOLYGON (((-73.89964 4...
8 8 3708 55881008 35002.58 MULTIPOLYGON (((-73.92478 4...
9 9 3709 124117798 73287.89 MULTIPOLYGON (((-73.83668 4...
10 10 3710 137760355 90067.74 MULTIPOLYGON (((-73.89681 4...
#for shapefiles, download the zipped shapefile directory into your current working directory, then unzip (can do this outside R if you prefer!)
#download
download.file(url = "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nypuma2010_21c.zip",destfile = "nypuma2010_21c.zip")
#unzip
unzip("nypuma2010_21c.zip")
#read in the shapefile
pumas_2010_shp <- st_read("nypuma2010_21c/nypuma2010.shp") %>% st_as_sf()
Reading layer `nypuma2010' from data source `/Users/sarahrankin/Documents/Mapping/nypuma2010_21c/nypuma2010.shp' using driver `ESRI Shapefile'
Simple feature collection with 55 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
Projected CRS: NAD83 / New York Long Island (ftUS)
pumas_2010_shp %>% print()
Simple feature collection with 55 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
Projected CRS: NAD83 / New York Long Island (ftUS)
First 10 features:
PUMA Shape_Leng Shape_Area geometry
1 3701 53227.11 97928517 MULTIPOLYGON (((1012885 268...
2 3702 106167.59 188993632 MULTIPOLYGON (((1021632 267...
3 3703 305269.14 267643637 MULTIPOLYGON (((1042822 243...
4 3704 47970.20 106216909 MULTIPOLYGON (((1026309 256...
5 3705 68697.60 122483690 MULTIPOLYGON (((1020081 255...
6 3706 51826.07 43886931 MULTIPOLYGON (((1018060 261...
7 3707 37374.60 42281072 MULTIPOLYGON (((1012009 253...
8 3708 35002.71 55881068 MULTIPOLYGON (((1005061 247...
9 3709 73288.78 124117768 MULTIPOLYGON (((1029456 237...
10 3710 90067.96 137760290 MULTIPOLYGON (((1012822 229...
geom_sf() to a ggplot renders the data in the geometry column - in this case it draws polygons because the geometry type is multipolygonpumas_2010_geojson %>%
ggplot() +
geom_sf()
geom_sf and geom_point calls ggplot() +
geom_sf(data = pumas_2010_geojson,
fill = "grey",
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency)) +
scale_color_manual(name = "System",values = library_system_colors) +
theme_void()
sf objects can have a CRS (coordinate reference system) defined
What if our point data doesn’t match our CRS?
#same plot but with the shapefile data rather than geoJSON
ggplot() +
geom_sf(data = pumas_2010_shp,
fill = "grey",
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency)) +
scale_color_manual(name = "System",values = library_system_colors) +
theme_void()
ggplot() +
geom_sf(data = pumas_2010_shp,
fill = "grey",
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=as.numeric(xcoord),
y=as.numeric(ycoord),
color = overagency)) +
scale_color_manual(name = "System",values = library_system_colors) +
theme_void()
st_transform and st_crs ggplot() +
geom_sf(data = pumas_2010_shp %>% st_transform(crs = st_crs(pumas_2010_geojson)),
fill = "grey",
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency)) +
scale_color_manual(name = "System",values = library_system_colors) +
theme_void()
NYC Open Data: NYC’s Internet Master Plan: Home Broadband and Mobile Broadband Adoption by PUMA
Based on American Community Survey Data; shows share of households that have both mobile and home broadband
broadband_use <- read_csv("https://data.cityofnewyork.us/api/views/g5ah-i2sh/rows.csv?accessType=DOWNLOAD")
── Column specification ────────────────────────────────────────────────────────────────────────────────────────
cols(
`PUMA (Public Use Microdata Sample Areas)` = col_double(),
Borough = col_character(),
`Home Broadband and Mobile Broadband Adoption (Percentage of Households)` = col_double(),
`Home Broadband and Mobile Broadband Adoption by Quartiles (High, Medium-High, Medium-Low, Low)` = col_character()
)
#clean up the names a bit
broadband_use <- broadband_use %>%
rename(PUMA = `PUMA (Public Use Microdata Sample Areas)`,
broadband_adoption = `Home Broadband and Mobile Broadband Adoption (Percentage of Households)`,
broadband_adoption_quartile = `Home Broadband and Mobile Broadband Adoption by Quartiles (High, Medium-High, Medium-Low, Low)`)
broadband_use %>% head()
#convert PUMA to string; make adoption quartile into factor so it'll sort correctly
broadband_use <- broadband_use %>%
mutate(PUMA = as.character(PUMA),
broadband_adoption_quartile = factor(broadband_adoption_quartile,
levels = c("High", "Medium High", "Medium Low", "Low")))
pumas_2010_geojson <- pumas_2010_geojson %>%
left_join(broadband_use, by = "PUMA")
pumas_2010_geojson %>% print()
Simple feature collection with 55 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91554
Geodetic CRS: WGS 84
First 10 features:
OBJECTID PUMA Shape__Area Shape__Length Borough broadband_adoption broadband_adoption_quartile
1 1 3701 97928415 53226.65 Bronx 0.46 High
2 2 3702 188993596 106167.48 Bronx 0.30 Low
3 3 3703 267645154 305260.19 Bronx 0.34 Medium Low
4 4 3704 106217077 47970.15 Bronx 0.30 Low
5 5 3705 122483734 68697.43 Bronx 0.32 Low
6 6 3706 43887042 51825.93 Bronx 0.39 Medium High
7 7 3707 42281072 37374.60 Bronx 0.42 High
8 8 3708 55881008 35002.58 Bronx 0.40 Medium High
9 9 3709 124117798 73287.89 Bronx 0.29 Low
10 10 3710 137760355 90067.74 Bronx 0.34 Medium Low
geometry
1 MULTIPOLYGON (((-73.89641 4...
2 MULTIPOLYGON (((-73.86477 4...
3 MULTIPOLYGON (((-73.78833 4...
4 MULTIPOLYGON (((-73.84793 4...
5 MULTIPOLYGON (((-73.87046 4...
6 MULTIPOLYGON (((-73.87773 4...
7 MULTIPOLYGON (((-73.89964 4...
8 MULTIPOLYGON (((-73.92478 4...
9 MULTIPOLYGON (((-73.83668 4...
10 MULTIPOLYGON (((-73.89681 4...
Map the geom_sf fill aesthetic to broadband_adoption
Use scale_fill_gradient to specify fill colors
ggplot() +
geom_sf(data = pumas_2010_geojson,
aes(fill = broadband_adoption),
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency),
size = 2) +
scale_color_manual(name = "System",values = library_system_colors) +
scale_fill_gradient(low = "grey70",high = "grey10", labels = scales::percent_format(accuracy = 1), name = "Broadband adoption") +
coord_sf() +
theme_void()
scale_fill_manualbroadband_cols <- c("grey10","grey30","grey50","grey70") %>% set_names(levels(pumas_2010_geojson$broadband_adoption_quartile))
ggplot() +
geom_sf(data = pumas_2010_geojson,
aes(fill = broadband_adoption_quartile),
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency),
size = 2) +
scale_color_manual(name = "System",values = library_system_colors) +
scale_fill_manual(values = broadband_cols, name = "Broadband adoption") +
coord_sf() +
theme_void()
https://www1.nyc.gov/site/planning/planning-level/nyc-population/american-community-survey.page
download.file(url = "https://www1.nyc.gov/assets/planning/download/office/planning-level/nyc-population/acs/econ_2018_acs5yr_puma.xlsx",
destfile = "econ_2018_acs5yr_puma.xlsx")
econdata_puma <- readxl::read_xlsx("econ_2018_acs5yr_puma.xlsx", sheet = "EconData")
pumas_2010_geojson <- pumas_2010_geojson %>%
left_join(econdata_puma %>%
select(PUMA = GeoID,GeogName,median_household_income = MdHHIncE) %>%
mutate(median_income_category =
cut(median_household_income,
breaks = c(0,50000,75000,100000,150000),
labels = c("$0-50K","$50-75K","$75-100K","$100K+"))))
Joining, by = "PUMA"
median_income_category_cols <- c("#a9ccbc","#7fb29b","#4c7f68","#335545") %>% set_names(levels(pumas_2010_geojson$median_income_category))
ggplot() +
geom_sf(data = pumas_2010_geojson,
aes(fill = median_income_category),
color = "grey"
) +
geom_point(data = public_libraries,
aes(x=longitude,
y=latitude,
color = overagency),
size = 2) +
scale_color_manual(name = "System",values = library_system_colors) +
scale_fill_manual(values = median_income_category_cols, name = "Median income") +
coord_sf() +
theme_void()
library(leaflet)
Leaflet: open source javascript library for making interactive maps on various platforms
R has a leaflet package that lets you create leaflet “widgets” within R
Create a map, add tiles, set the view
Similar to ggplot in that we keep adding layers/elements, except we continue using the magrittr pipe %>% to add elements to a leaflet map (instead of the ggplot +)
for now we’re creating a map without any data
leaflet() %>%
#addTiles() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10)
Add our public library branches
addCircles() adds points to the map
can specify the data source within the addCircles() call or in the original leaflet() call
to refer to fields in our data frame, use the tilde prefix ~
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10) %>%
addCircles(data = public_libraries,
lng = ~longitude,
lat = ~latitude
)
Use colorFactor() to create a palette to color the circles by system
Use the label and popup arguments to add hover/click info
pal_system <- colorFactor(palette = library_system_colors ,
domain = unique(public_libraries$overagency),
ordered = T)
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10) %>%
addCircles(data = public_libraries,
lng = ~longitude,
lat = ~latitude,
color = ~pal_system(overagency),
fill = ~pal_system(overagency),
radius = 100,
label = ~facname,
popup = ~address
) %>%
addLegend(position = "topleft",
pal = pal_system,
values = unique(public_libraries$overagency))
Because we have a basemap, we don’t need the PUMA polygons to create a basic map. But we can add them if we want our chloropleths on the interactive map.
Fill color argument in leaflet is fillColor; specify opacity with fillOpacity
Use stroke to specify whether boundary lines are drawn; color , weight, and opacity to control their appearance
pal_broadband <- colorFactor(palette = broadband_cols,
domain = levels(pumas_2010_geojson$broadband_adoption_quartile),
ordered = T)
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10) %>%
addPolygons(data = pumas_2010_geojson,
fillColor = ~pal_broadband(broadband_adoption_quartile),
fillOpacity = .7,
color = "grey",
stroke = T,
weight = 1,
label = ~PUMA,
popup = ~paste0("<b>PUMA ",PUMA,"</b>",
"<br>", GeogName,
"<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
"<br>Broadband adoption category: ",broadband_adoption_quartile))
Combine it with the libraries and add a legend
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10) %>%
addPolygons(data = pumas_2010_geojson,
fillColor = ~pal_broadband(broadband_adoption_quartile),
fillOpacity = .7,
color = "grey",
stroke = T,
weight = 1,
label = ~PUMA,
popup = ~paste0("<b>PUMA ",PUMA,"</b>",
"<br>", GeogName,
"<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
"<br>Broadband adoption category: ",broadband_adoption_quartile)) %>%
addCircles(data = public_libraries,
lng = ~longitude,
lat = ~latitude,
color = ~pal_system(overagency),
fill = ~pal_system(overagency),
radius = 100,
label = ~facname,
popup = ~address
) %>%
addLegend(position = "topleft",
pal = pal_system,
values = public_libraries$overagency,
title = "Library System") %>%
addLegend(position = "topleft",
pal = pal_broadband,
values = pumas_2010_geojson$broadband_adoption_quartile,
opacity = .7,
title = "Broadband Access")
Visualize more than one layer
Use group argument within each layer to specify which group it belongs to
addLayersControl adds a control to toggle between layers or turn layers on and off
assigning the map to an object creates a leaflet object with class htmlwidget
pal_income <- colorFactor(palette = median_income_category_cols,
domain = levels(pumas_2010_geojson$median_income_category),
ordered = T)
libraries_map <- leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
setView(lat = 40.7,lng = -74, zoom = 10) %>%
addPolygons(data = pumas_2010_geojson,
fillColor = ~pal_broadband(broadband_adoption_quartile),
fillOpacity = .7,
color = "grey",
stroke = T,
weight = 1,
label = ~PUMA,
popup = ~paste0("<b>PUMA ",PUMA,"</b>",
"<br>", GeogName,
"<br>Broadband adoption: ",round(broadband_adoption*100,0),"%",
"<br>Broadband adoption category: ",broadband_adoption_quartile),
group = "Broadband") %>%
addPolygons(data = pumas_2010_geojson,
fillColor = ~pal_income(median_income_category),
fillOpacity = .7,
color = "grey",
stroke = T,
weight = 1,
label = ~PUMA,
popup = ~paste0("<b>PUMA ",PUMA,"</b>",
"<br>", GeogName,
"<br>Median income: $",format(median_household_income,big.mark = ","),
"<br>Median income category: ",median_income_category),
group = "income") %>%
addCircles(data = public_libraries,
lng = ~longitude,
lat = ~latitude,
color = ~pal_system(overagency),
fill = ~pal_system(overagency),
opacity = 1,
fillOpacity = 1,
radius = 100,
label = ~facname,
popup = ~address#,
#group = "Broadband"
) %>%
# addCircles(data = public_libraries,
# lng = ~longitude,
# lat = ~latitude,
# color = ~pal_system(overagency),
# fill = ~pal_system(overagency),
# opacity = 1,
# fillOpacity = 1,
# radius = 100,
# label = ~facname,
# popup = ~address,
# group = "Income"
# ) %>%
addLegend(position = "topleft",
pal = pal_system,
values = public_libraries$overagency,
title = "Library System") %>%
addLegend(position = "topleft",
pal = pal_broadband,
values = pumas_2010_geojson$broadband_adoption_quartile,
opacity = .7,
title = "Broadband Access",
group = "Broadband") %>%
addLegend(position = "topleft",
pal = pal_income,
values = pumas_2010_geojson$median_income_category,
opacity = .7,
title = "Median household income",
group = "Income") %>%
addLayersControl(
baseGroups = c("Broadband","Income"),
options = layersControlOptions(collapsed = F)#,
#overlayGroups = c("Branches")
)
libraries_map
The htmlwidgets package lets you save this map as an html file, which can then be opened in a browser or embedded in a web page
htmlwidgets::saveWidget(libraries_map,
file = "libraries_map.html",
title = "NYC Libraries: Broadband Access and Median Income")
tidycensus package: amazing R interface for Census/ACS API, including geometries
tigris package: more Census geometries
rnaturalearth: world map data
tmap: Alternative package for tmapping in R, with ggplot-like syntax - some find this the best way to get started
mapboxapi package: more basemap options, geocoding services, isochrones (requires setting up a mapbox account)
Geocoding
Within NYC, Geosupport